Summary: calibrating primer sets from a real experimental test
This vignette shows how to use tidyqpcr functions to calibrate qPCR probes.
This is real qPCR data by WD and XH in July 2023, testing new RT-qPCR primer sets against S. cerevisiae genes. We took exponential-phase total RNA previously extracted by Jamie Auxillos.
We tested 3 primer sets each for 3 genes: * NOB1/YOR056C (3 primer sets) * FCP1/YMR277C (3 primer sets) * ARG8/YOL140W (3 primer sets) * PGK1 (1 primer set as positivel control)
We started with two biological replicate RNA samples, treated with DNase and then split for a test sample with reverse transcriptase (RT) and negative control without reverse transcriptase (-RT). We also took a no template (NT) negative control. For each RT reaction we do serial 5x dilutions down to 125x to form a quantitative calibration curve.
The data were measured on a Roche LC480 instrument in a single 384-well plate. Quantification was performed in the Roche LightCycler software prior to exporting the data.
Setup knitr options and load packages
# knitr options for report generation
knitr::opts_chunk$set(
warning = FALSE, message = FALSE, echo = TRUE, cache = FALSE,
results = "show"
)
# Load packages
library(tidyr)
library(ggplot2)
library(dplyr)
library(tidyqpcr)
# set default theme for graphics
theme_set(theme_bw(base_size = 11) %+replace%
theme(
strip.background = element_blank(),
panel.grid = element_blank()
))Set up experiment
Describe which primer set we put in which well using a row key
In this experiment, each primer set was in a different row of the 384-well plate. We describe this by creating a row key, a data frame describing the rows of the plate, the primer sets, and the genes that they detect.
We refer to each primer set as a target_id, because each
primer set has a different target amplicon, on a different location in a
gene. tidyqpcr insists on having a variable called
target_id that uniquely identifies each different target
that you detect.
# Names of target genes
gene_name_levels <- c("NOB1", "FCP1", "ARG8", "PGK1")
# ORF ids of target genes
target_levels <- c("YOR056C", "YMR277C", "YOL140W","positive_control")
# Repeats of gene names to account for testing multiple primer sets
gene_name_values <- c(rep(gene_name_levels[1:3], each = 3),
rep(gene_name_levels[4], each = 1))
# id numbers of multiple probesets (reflecting IDs as ordered)
target_id_levels <- paste(gene_name_values,
c(1, 2, 3, 1, 2, 3, 1, 2, 3, 1),
sep = "_"
)
# 4 is positive control primer
rowkey <- tibble(
well_row = LETTERS[1:10],
gene_name = factor(gene_name_values, levels = gene_name_levels),
target_id = factor(target_id_levels, levels = target_id_levels)
)
print(rowkey)## # A tibble: 10 × 3
## well_row gene_name target_id
## <chr> <fct> <fct>
## 1 A NOB1 NOB1_1
## 2 B NOB1 NOB1_2
## 3 C NOB1 NOB1_3
## 4 D FCP1 FCP1_1
## 5 E FCP1 FCP1_2
## 6 F FCP1 FCP1_3
## 7 G ARG8 ARG8_1
## 8 H ARG8 ARG8_2
## 9 I ARG8 ARG8_3
## 10 J PGK1 PGK1_1Combine the row key describing primer sets with column key describing on samples and dilutions.
We use a default design built in to tidyqpcr,
create_colkey_4diln_2ctrl_in_24.
plate1plan <-
label_plate_rowcol(
create_blank_plate(),
rowkey,
create_colkey_4diln_2ctrl_in_24()
) %>%
mutate(sample_id = paste(biol_rep, dilution_nice, sep = "_"))Spot-check the plate plan
Checks that for selected technical replicate/probe/dilution combinations, the plate contains the right number of replicates.
plate1plan %>%
filter(tech_rep == "1",
target_id == target_id_levels[1],
dilution_nice == "1x")## # A tibble: 2 × 11
## well well_row well_col dilution dilution_nice prep_type biol_rep tech_rep
## <chr> <fct> <fct> <dbl> <chr> <fct> <fct> <fct>
## 1 A1 A 1 1 1x +RT A 1
## 2 A13 A 13 1 1x +RT B 1
## # ℹ 3 more variables: gene_name <fct>, target_id <fct>, sample_id <chr>plate1plan %>%
filter(tech_rep == "2",
target_id == target_id_levels[4])## # A tibble: 12 × 11
## well well_row well_col dilution dilution_nice prep_type biol_rep tech_rep
## <chr> <fct> <fct> <dbl> <chr> <fct> <fct> <fct>
## 1 D7 D 7 1 1x +RT A 2
## 2 D8 D 8 0.2 5x +RT A 2
## 3 D9 D 9 0.04 25x +RT A 2
## 4 D10 D 10 0.008 125x +RT A 2
## 5 D11 D 11 1 -RT -RT A 2
## 6 D12 D 12 1 NT NT A 2
## 7 D19 D 19 1 1x +RT B 2
## 8 D20 D 20 0.2 5x +RT B 2
## 9 D21 D 21 0.04 25x +RT B 2
## 10 D22 D 22 0.008 125x +RT B 2
## 11 D23 D 23 1 -RT -RT B 2
## 12 D24 D 24 1 NT NT B 2
## # ℹ 3 more variables: gene_name <fct>, target_id <fct>, sample_id <chr>Display the plate plan
This can be printed out to facilitate loading the plate correctly.
display_plate_qpcr(plate1plan)
Analyse Cq (quantification cycle count) data
Load and summarize data
# NOTE: system.file() accesses data from this R package
# To use your own data, remove the call to system.file(),
# instead pass your data's filename to read_lightcycler_1colour_cq()
# or to another relevant read_ function
plates <-
read_lightcycler_1colour_cq("WD_05072023_primers_ct.txt") %>%
right_join(plate1plan, by = "well")plates## # A tibble: 384 × 18
## include color well sample_info cq concentration standard status well_row
## <lgl> <int> <chr> <chr> <dbl> <dbl> <int> <lgl> <fct>
## 1 TRUE 255 A1 Sample 1 16.4 NA 0 NA A
## 2 TRUE 255 A2 Sample 2 18.7 NA 0 NA A
## 3 TRUE 255 A3 Sample 3 21.2 NA 0 NA A
## 4 TRUE 255 A4 Sample 4 23.7 NA 0 NA A
## 5 TRUE 255 A5 Sample 5 33.7 NA 0 NA A
## 6 TRUE 255 A6 Sample 6 38.6 NA 0 NA A
## 7 TRUE 255 A7 Sample 7 16.1 NA 0 NA A
## 8 TRUE 255 A8 Sample 8 18.3 NA 0 NA A
## 9 TRUE 255 A9 Sample 9 20.7 NA 0 NA A
## 10 TRUE 255 A10 Sample 10 23.2 NA 0 NA A
## # ℹ 374 more rows
## # ℹ 9 more variables: well_col <fct>, dilution <dbl>, dilution_nice <chr>,
## # prep_type <fct>, biol_rep <fct>, tech_rep <fct>, gene_name <fct>,
## # target_id <fct>, sample_id <chr>summary(plates)## include color well sample_info
## Mode:logical Min. : 255 Length:384 Length:384
## TRUE:384 1st Qu.: 255 Class :character Class :character
## Median : 255 Mode :character Mode :character
## Mean :32090
## 3rd Qu.:65280
## Max. :65280
##
## cq concentration standard status well_row
## Min. : 9.15 Min. : NA Min. :0 Mode:logical A : 24
## 1st Qu.:17.22 1st Qu.: NA 1st Qu.:0 NA's:384 B : 24
## Median :20.40 Median : NA Median :0 C : 24
## Mean :21.65 Mean :NaN Mean :0 D : 24
## 3rd Qu.:24.46 3rd Qu.: NA 3rd Qu.:0 E : 24
## Max. :38.59 Max. : NA Max. :0 F : 24
## NA's :188 NA's :384 (Other):240
## well_col dilution dilution_nice prep_type biol_rep tech_rep
## 1 : 16 Min. :0.0080 Length:384 +RT:256 A:192 1:192
## 2 : 16 1st Qu.:0.0400 Class :character -RT: 64 B:192 2:192
## 3 : 16 Median :0.6000 Mode :character NT : 64
## 4 : 16 Mean :0.5413
## 5 : 16 3rd Qu.:1.0000
## 6 : 16 Max. :1.0000
## (Other):288
## gene_name target_id sample_id
## NOB1: 72 NOB1_1 : 24 Length:384
## FCP1: 72 NOB1_2 : 24 Class :character
## ARG8: 72 NOB1_3 : 24 Mode :character
## PGK1: 24 FCP1_1 : 24
## NA's:144 FCP1_2 : 24
## (Other):120
## NA's :144Visualise Cq values for each well.
Visualising the Cq values shows that the Cq value is different for
each primer set in each row. Within each section of a row for a single
replicate of dilutions, Cq consistently increases with dilutions as
expected. The grey tiles for most -RT and NT columns mean that the value
is NA, i.e. no Cq value was reported. This is good.
display_plate_value(plates)
Visualisation might also help to identify unwanted positional effects. For example, a PCR machine is broken, wells close to an edge of the plate can have different behaviour from central wells.
Plot unnormalized data shows that -RT and NT controls are low
This plot visualises the Cq data in a way that highlights the meaning instead of the position on the plate. Again, it shows that the Cq value is different for each primer set, and that for each primer st Cq consistently increases with dilutions as expected.
Again, we detect no signal in NT (no template) negative control so those points are mostly missing. There is a very weak signal with high Cq in some -RT (no reverse transcriptase) negative controls.
ggplot(data = plates) +
geom_point(aes(x = target_id,
y = cq,
colour = dilution_nice,
shape = prep_type),
position = position_jitter(width = 0.2, height = 0)
) +
labs(
y = "Cycle count to threshold",
title = "All reps, unnormalized"
) +
facet_wrap(~biol_rep) +
scale_y_continuous(breaks = seq(from = 15, to = 35, by = 5)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
panel.grid.major.y = element_line(colour="grey80", size = 0.2))
Dilution series is linear for all probes
Visual display of linearity of cq with log(dilution).
ggplot(data = filter(plates, prep_type == "+RT"), aes(x = dilution, y = cq)) +
geom_point() +
stat_smooth(
formula = y ~ x, method = "lm", se = FALSE,
aes(colour = "fit", linetype = "fit")
) +
stat_smooth(
formula = y ~ 1 + offset(-x * log(10) / log(2)), method = "lm", se = FALSE,
aes(colour = "theory", linetype = "theory")
) +
scale_x_log10(breaks = 10 ^ - (0:3)) +
scale_y_continuous(breaks = seq(0, 30, 2)) +
labs(
y = "Cycle count to threshold",
title = "All reps, unnormalized",
colour = "Dilution", linetype = "Dilution"
) +
facet_grid(target_id ~ biol_rep, scales = "free_y") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5))
Calculate primer efficiencies for all probes
Use regression to estimate linearity of cq with log(dilution), including the slope or efficiency.
# calculate_efficiency_bytargetid(plates)
# failedDilution series for nice probes only shows linearity clearly
target_id_levels_niceprobes <-
c("NOB1_3", "FCP1_2", "FCP1_3", "ARG8_1",
"ARG8_3")
ggplot(
data = filter(plates,
prep_type == "+RT",
target_id %in% target_id_levels_niceprobes),
aes(x = dilution, y = cq)
) +
geom_point() +
stat_smooth(
formula = y ~ x, method = "lm", se = FALSE,
aes(colour = "fit", linetype = "fit")
) +
stat_smooth(
formula = y ~ 1 + offset(-x * log(10) / log(2)),
method = "lm",
se = FALSE,
aes(colour = "theory", linetype = "theory")
) +
scale_x_log10(breaks = 10 ^ - (0:3)) +
scale_y_continuous(breaks = seq(0, 30, 2)) +
labs(
y = "Cycle count to threshold",
title = "All reps, unnormalized",
colour = "Dilution", linetype = "Dilution"
) +
facet_grid(target_id ~ biol_rep, scales = "free_y") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5))
Analyse amplification and melt curve data
Load raw data for amplification and melt curves.
# NOTE: system.file() accesses data from this R package
# To use your own data, remove the call to system.file(),
# instead pass your data's filename to read_lightcycler_1colour_cq()
# or to another relevant read_ function
plate1curve <-
read_lightcycler_1colour_raw("WD_05072023_primers.txt") %>%
debaseline() %>%
left_join(plate1plan, by = "well")
# amplification curve is program 2
platesamp <- plate1curve %>%
filter(program_no == 2)
# melt curve is program 3 or 4, depending on cycler setting
platesmelt <- plate1curve %>%
filter(program_no == 3) %>%
calculate_drdt_plate() %>%
filter(temperature >= 61)Plot de-baseline’d raw data for single well
ggplot(
data = platesamp %>% filter(well == "A1"),
aes(x = cycle, y = fluor_signal)
) +
geom_line() +
scale_y_continuous(expand = c(0.01, 0.01))
Plot all amplification curves
Broken up by technical replicate here, to avoid overplotting.
ggplot(
data = platesamp %>%
filter(tech_rep == "1"),
aes(x = cycle,
y = fluor_signal,
colour = factor(dilution),
linetype = prep_type)
) +
facet_grid(target_id ~ biol_rep, scales = "free_y") +
scale_linetype_manual(values = c("+RT" = "solid",
"-RT" = "dashed",
"NT" = "dotted")) +
geom_line() +
scale_x_continuous(breaks = seq(60, 100, 10),
minor_breaks = seq(60, 100, 5)) +
labs(title = "All Amp Curves, tech_rep A")
ggplot(
data = platesamp %>%
filter(tech_rep == "2"),
aes(x = cycle,
y = fluor_signal,
colour = factor(dilution),
linetype = prep_type)
) +
facet_grid(target_id ~ biol_rep, scales = "free_y") +
scale_linetype_manual(values = c("+RT" = "solid",
"-RT" = "dashed",
"NT" = "dotted")) +
geom_line() +
scale_x_continuous(breaks = seq(60, 100, 10),
minor_breaks = seq(60, 100, 5)) +
labs(title = "All Amp Curves, tech_rep B")
Plot melt curve for single well
ggplot(
data = platesmelt %>%
filter(well == "A1"),
aes(x = temperature, y = dRdT)
) +
facet_wrap(~target_id) +
geom_line() +
scale_y_continuous(expand = c(0.02, 0.02))
Plot all melt curves
Again broken up by technical replicate.
ggplot(
data = platesmelt %>%
filter(tech_rep == "1"),
aes(x = temperature,
y = dRdT,
colour = factor(dilution),
linetype = prep_type)
) +
facet_grid(target_id ~ biol_rep, scales = "free_y") +
scale_linetype_manual(values = c("+RT" = "solid",
"-RT" = "dashed",
"NT" = "dotted")) +
geom_line() +
scale_x_continuous(breaks = seq(60, 100, 10),
minor_breaks = seq(60, 100, 5)) +
labs(title = "All Melt Curves, tech_rep A")
ggplot(
data = platesmelt %>%
filter(tech_rep == "2"),
aes(x = temperature,
y = dRdT,
colour = factor(dilution),
linetype = prep_type)
) +
facet_grid(target_id ~ biol_rep, scales = "free_y") +
scale_linetype_manual(values = c("+RT" = "solid",
"-RT" = "dashed",
"NT" = "dotted")) +
geom_line() +
scale_x_continuous(breaks = seq(60, 100, 10),
minor_breaks = seq(60, 100, 5)) +
labs(title = "All Melt Curves, tech_rep B")
Plot zoomed melt curves
ggplot(
data = platesmelt %>%
filter(tech_rep == "1", prep_type == "+RT"),
aes(x = temperature, y = dRdT, colour = factor(dilution))
) +
facet_grid(target_id ~ biol_rep, scales = "free_y") +
geom_line() +
scale_x_continuous(
breaks = seq(60, 100, 5),
minor_breaks = seq(60, 100, 1),
limits = c(73, 87)
) +
labs(title = "Melt curves, zoomed, tech_rep A") +
theme(
panel.grid.major.x = element_line(colour = "grey50", size = 0.4),
panel.grid.minor.x = element_line(colour = "grey70", size = 0.1)
)
ggplot(
data = platesmelt %>%
filter(tech_rep == "2", prep_type == "+RT"),
aes(x = temperature, y = dRdT, colour = factor(dilution))
) +
facet_grid(target_id ~ biol_rep, scales = "free_y") +
geom_line() +
scale_x_continuous(
breaks = seq(60, 100, 5),
minor_breaks = seq(60, 100, 1),
limits = c(73, 87)
) +
labs(title = "Melt curves, zoomed, tech_rep B") +
theme(
panel.grid.major.x = element_line(colour = "grey50", size = 0.4),
panel.grid.minor.x = element_line(colour = "grey70", size = 0.1)
)
Plot only zoomed melt curves for nice probes
ggplot(
data = platesmelt %>%
filter(
tech_rep == "1",
prep_type == "+RT",
dilution_nice == "1x",
target_id %in% target_id_levels_niceprobes
),
aes(x = temperature, y = dRdT, colour = biol_rep)
) +
facet_grid(target_id ~ ., scales = "free_y") +
geom_line() +
scale_x_continuous(
breaks = seq(60, 100, 5),
minor_breaks = seq(60, 100, 1),
limits = c(73, 87)
) +
labs(title = "Nice probes, tech_rep A") +
theme(
panel.grid.major.x = element_line(colour = "grey50", size = 0.4),
panel.grid.minor.x = element_line(colour = "grey70", size = 0.1)
)
Conclude acceptable primer sets
These probes have sensible standard curves, amplification curves, melt curves. In tie-breakers we pick the more highly detected probe.